Linear Mixed Models (LMMs) - Moderation

Joshua F. Wiley

2020-05-31

Download the raw R markdown code here https://jwiley.github.io/MonashHonoursStatistics/LMM_Moderation.rmd. These are the R packages we will use.

options(digits = 4)

## emmeans is a new package

library(data.table)
## data.table 1.12.8 using 24 threads (see ?getDTthreads).  Latest news: r-datatable.com
library(JWileymisc)
library(extraoperators)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(multilevelTools)
library(visreg)
library(ggplot2)
library(ggpubr)
## Loading required package: magrittr
library(haven)
library(emmeans)

## load data collection exercise data
## merged is a a merged long dataset of baseline and daily
dm <- as.data.table(read_sav("Merged.sav"))

1 LMM Notation

Let’s consider the formula for a relatively simple LMM:

\[ y_{ij} = b_{0j} + b_1 * x_{1j} + b_2 * x_{2ij} + \varepsilon_{ij} \]

Here as before, the i indicates the ith observation for a specific unit (e.g., person but the unit could also be classrooms, doctors, etc.) and the j indicates the jth unit (in psychology usually person).

Regression coefficients, the \(b\)s with a j subscript indicate a fixed and random effect. That is, the coefficient is allowed to vary across the units, j. As before, these coefficients in practice are decomposed into a fixed and random part:

\[ b_{0j} = \gamma_{00} + u_{0j} \]

and we estimate in our LMM the fixed effect part, \(\gamma_{00}\), and the variance / standard deviation of the random effect or the covariance matrix if there are multiple random effects, \(\mathbf{G}\):

\[ u_{0j} \sim \mathcal{N}(0, \mathbf{G}) \]

Regression coefficients without any j subscript indicate fixed only effects, effects that do not vary across units, j. These are fixed effects and get estimated directly.

Predictors / explanatory variables, the \(x\)s with an i subscript indicate that the variable varies within a unit. Note that the outcome, \(y\) must vary within units to be used in a LMM.

In this case, the notation tells us the following:

The following decision tree provides some guide to when a predictor / explanatory variable can be a fixed and random effect.

Type of effect decision tree

Let’s see two examples of putting this basic model into practice.

\[ mood_{ij} = b_{0j} + b_1 * age_{j} + b_2 * stress_{ij} + \varepsilon_{ij} \]

The corresponding R code is:

summary(lmer(mood ~ age + stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 573.9
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.377 -0.464  0.188  0.587  2.301 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.455    0.675   
##  Residual             1.003    1.001   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   6.5149     1.0067  62.4720    6.47  1.7e-08 ***
## age          -0.0217     0.0415  53.9863   -0.52      0.6    
## stress       -0.2838     0.0652 176.8111   -4.35  2.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) age   
## age    -0.962       
## stress -0.329  0.088

Here is another example decomposing stress into a between and within component.

\[ mood_{ij} = b_{0j} + b_1 * Bstress_{j} + b_2 * Wstress_{ij} + \varepsilon_{ij} \]

dm[, c("Bstress", "Wstress") := meanDeviations(stress), by = ID]

summary(lmer(energy ~ Bstress + Wstress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: energy ~ Bstress + Wstress + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 595.4
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.763 -0.580  0.112  0.702  1.835 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.423    0.65    
##  Residual             1.187    1.09    
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   5.3788     0.4622  48.0445   11.64  1.4e-15 ***
## Bstress      -0.3008     0.1170  47.1779   -2.57    0.013 *  
## Wstress      -0.2100     0.0855 139.0037   -2.46    0.015 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) Bstrss
## Bstress -0.962       
## Wstress  0.000  0.000

We can make more effects random effects. For example, taking our earlier example and just changing \(b_2\) into \(b_{2j}\):

\[ mood_{ij} = b_{0j} + b_1 * age_{j} + b_{2j} * stress_{ij} + \varepsilon_{ij} \]

The corresponding R code is:

summary(lmer(mood ~ age + stress + (stress | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + (stress | ID)
##    Data: dm
## 
## REML criterion at convergence: 570.7
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.525 -0.510  0.188  0.588  2.361 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 0.7759   0.881         
##           stress      0.0644   0.254    -0.78
##  Residual             0.9388   0.969         
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   6.3282     0.9345 43.1574    6.77  2.7e-08 ***
## age          -0.0131     0.0380 36.7556   -0.34   0.7329    
## stress       -0.2762     0.0758 27.3002   -3.64   0.0011 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) age   
## age    -0.951       
## stress -0.346  0.067

technically note now with two random effects, we assume that the random effects, \(u_{0j}\) and \(u_{2j}\), which we collectively denote \(\mathbf{u}_{j}\) follow a multivariate normal distribution with covariance matrix \(\mathbf{G}\).

\[ \mathbf{u}_{j} \sim \mathcal{N}(0, \mathbf{G}) \]

Based on the little decision chart, between unit only variables, like \(age_j\) and \(Bstress_j\) cannot be random effects. Also, while it is technically possible for something to only be a random effect without a corresponding fixed effect, its not common and not recommended as it would be equivalent to assuming that the fixed effect, the mean of the distribution, is 0, which is rarely appropriate.

2 Interactions in LMMs

Interactions in LMMs work effectively the same way that interactions in GLMs do, although there are a few nuances in options and possible interpretations. Using the notation from above, let’s consider a few different possible interactions.

2.1 Cross Level (Between and Within Unit) Interactions

First, let’s take our model with age and stress and include an interaction. Here is the model without an interaction.

\[ mood_{ij} = b_{0j} + b_1 * age_{j} + b_2 * stress_{ij} + \varepsilon_{ij} \]

The corresponding R code is:

summary(lmer(mood ~ age + stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 573.9
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.377 -0.464  0.188  0.587  2.301 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.455    0.675   
##  Residual             1.003    1.001   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   6.5149     1.0067  62.4720    6.47  1.7e-08 ***
## age          -0.0217     0.0415  53.9863   -0.52      0.6    
## stress       -0.2838     0.0652 176.8111   -4.35  2.3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) age   
## age    -0.962       
## stress -0.329  0.088

Now let’s add the interaction, as a fixed effect.

\[ mood_{ij} = b_{0j} + b_1 * age_{j} + b_2 * stress_{ij} + b_3 * (age_{j} * stress_{ij}) + \varepsilon_{ij} \]

The corresponding R code is:

## long way
summary(lmer(mood ~ age + stress + age:stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + age:stress + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 576.1
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.374 -0.443  0.163  0.580  2.311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.433    0.658   
##  Residual             0.996    0.998   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)   2.6818     2.3243 176.9101    1.15     0.25  
## age           0.1476     0.1014 177.7585    1.46     0.15  
## stress        0.7967     0.5968 168.1832    1.33     0.18  
## age:stress   -0.0480     0.0263 165.4710   -1.82     0.07 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) age    stress
## age        -0.993              
## stress     -0.914  0.914       
## age:stress  0.904 -0.915 -0.994
## short hand in R for simple main effect + interaction
## identical, but shorter to the above
summary(lmer(mood ~ age * stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age * stress + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 576.1
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.374 -0.443  0.163  0.580  2.311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.433    0.658   
##  Residual             0.996    0.998   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)   2.6818     2.3243 176.9101    1.15     0.25  
## age           0.1476     0.1014 177.7585    1.46     0.15  
## stress        0.7967     0.5968 168.1832    1.33     0.18  
## age:stress   -0.0480     0.0263 165.4710   -1.82     0.07 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) age    stress
## age        -0.993              
## stress     -0.914  0.914       
## age:stress  0.904 -0.915 -0.994

The relevant, new, part is the interaction term, \(b_3\), a fixed effect in this case. If we focus just on that one term, we see that the coefficient, \(b_3\) is applied to the arithmetic product of two variables, here age and stress. As it happens, one of them, age, varies only between units whereas the other, stress, varies within units. You will sometimes see this termed as “cross level” interaction between it involves a between and within varying variable.

\[ b_3 * (age_{j} * stress_{ij}) \]

As with interactions for regular GLMs, interactions in LMMs can be interpretted in different ways. The two common interpretations are easiest to see by factoring the regression equation. Here are three equal equations that highlight different ways of viewing the interaction.

In the latter two formats, it highlights how the simple effect of stress varies by age and how the simple effect of age varies by stress.

\[ \begin{align} mood_{ij} &= b_{0j} + b_1 * age_{j} + b_2 * stress_{ij} + b_3 * (age_{j} * stress_{ij}) + \varepsilon_{ij} \\ &= b_{0j} + b_1 * age_{j} + (b_2 + b_3 * age_j) * stress_{ij} + \varepsilon_{ij} \\ &= b_{0j} + (b_1 + b_3 * stress_{ij}) * age_{j} + b_2 * stress_{ij} + \varepsilon_{ij} \\ \end{align} \]

The nuance in LMMs comes in because some variables vary only between units and others within units. For example, when interpretting the interaction with respect to the simple effect of stress, we could say that the association between daily stress and mood on the same day depends on the age of a participant. Conversely, when interpretting with respect to the simple effect of age, we could say that the association of participant age and average mood depends on how stressed someone is feeling on a given day. Age varies only between people, stress varies within people, so that must be taken into account in the interpretation.

2.2 Between Unit Interactions

The same approach would work with other type of variables in LMMs. For example, here we have a model with age and female as predictors. Both vary only between units.

\[ \begin{align} mood_{ij} &= b_{0j} + b_1 * age_{j} + b_2 * female_{j} + b_3 * (age_{j} * female_{j}) + \varepsilon_{ij} \\ &= b_{0j} + b_1 * age_{j} + (b_2 + b_3 * age_j) * female_{j} + \varepsilon_{ij} \\ &= b_{0j} + (b_1 + b_3 * female_{j}) * age_{j} + b_2 * female_{j} + \varepsilon_{ij} \\ \end{align} \]

When interpretting the interaction with respect to the simple effect of female, we could say that the association between participant sex and average mood depends on the age of a participant. Conversely, when interpretting with respect to the simple effect of age, we could say that the association of participant age and average mood depends on participant sex.

2.3 Within Unit Interactions

Finally, both variables could vary within units.

\[ \begin{align} mood_{ij} &= b_{0j} + b_1 * energy_{ij} + b_2 * stress_{ij} + b_3 * (energy_{ij} * stress_{ij}) + \varepsilon_{ij} \\ &= b_{0j} + b_1 * energy_{ij} + (b_2 + b_3 * energy_{ij}) * stress_{ij} + \varepsilon_{ij} \\ &= b_{0j} + (b_1 + b_3 * stress_{ij}) * energy_{ij} + b_2 * stress_{ij} + \varepsilon_{ij} \\ \end{align} \]

When interpretting the interaction with respect to the simple effect of stress, we could say that the association between daily stress and mood on the same day depends on same day energy level. Conversely, when interpretting with respect to the simple effect of energy, we could say that the association of daily energy and same day mood depends on how stressed someone is feeling on a given day.

2.4 ‘Prediction’ Interpretation of Interactions

When one variable in an interaction varies within units and particularly when it is a random effect, another way that people sometimes interpret interactions is that the other variable is ‘predicting’ the random effect. This occurs most often when the moderator variable varies between units only. An example may be clearer than words.

First, we fit a model with age and stress as fixed effects predictors and a random slope as well for stress, stored in m1. Then we fit the same model but adding a fixed effect interaction between stress (within) and age (between). Calculating the difference in the variance of the random stress slope yields a sort of \(R^2\) measure of the variance in the random slope explained by the stress x age interaction.

The corresponding R code is:

## main effects only
m1 <- lmer(mood ~ age + stress + (1 + stress | ID), data = dm, REML=FALSE)

## interaction model
m2 <- lmer(mood ~ age * stress + (1 + stress | ID), data = dm, REML=FALSE)

## summary of both models, get random slope for stress variance
summary(m1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mood ~ age + stress + (1 + stress | ID)
##    Data: dm
## 
##      AIC      BIC   logLik deviance df.resid 
##    574.2    596.6   -280.1    560.2      176 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.581 -0.511  0.215  0.589  2.356 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 0.6643   0.815         
##           stress      0.0587   0.242    -0.77
##  Residual             0.9388   0.969         
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   6.2889     0.9103 44.2464    6.91  1.5e-08 ***
## age          -0.0115     0.0370 37.5843   -0.31  0.75696    
## stress       -0.2754     0.0741 27.4601   -3.72  0.00092 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) age   
## age    -0.951       
## stress -0.346  0.068
summary(m2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mood ~ age * stress + (1 + stress | ID)
##    Data: dm
## 
##      AIC      BIC   logLik deviance df.resid 
##    572.8    598.5   -278.4    556.8      175 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.628 -0.375  0.201  0.596  2.365 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID       (Intercept) 0.5051   0.711         
##           stress      0.0497   0.223    -0.73
##  Residual             0.9361   0.968         
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)   2.3520     2.2723 44.8968    1.04    0.306  
## age           0.1620     0.0993 45.3903    1.63    0.110  
## stress        0.8907     0.6317 53.3680    1.41    0.164  
## age:stress   -0.0517     0.0279 54.2982   -1.85    0.069 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) age    stress
## age        -0.993              
## stress     -0.929  0.928       
## age:stress  0.920 -0.931 -0.994
## variance in random stress slope explained by age
(.0587 - .0497) / .0587
## [1] 0.1533

Because age and stress interact in the model (as a fixed effect), the average stress slope no longer has to be the same for everyone. It can differ depending on their age. Specifically, the predicted average (fixed effects) stress slope for the jth person is \(b_{2j} + b_3 * age_j\). Recall that we normally break up random slopes into a fixed and random part:

\[ b_{2j} = \gamma_{20} + u_{2j} \]

Without the interaction by age, any participant level differences must go into the \(u_{2j}\) component, on which the variance/standard deviation of the slope is calculated. Now the simple stress slope for a given participant is:

\[ \gamma_{20} + b_3 * age_j + u_{2j} \]

so now the \(u_{2j}\) component, on which the variance/standard deviation of the slope is calculated, captures deviations from the fixed part, which includes both the average stress slope and a modification based on participant age. To the extent that \(b_3\) is different from 0, this will essentially reduce some differences that otherwise go into the \(u_{2j}\) component and thus will explain some of the variance in the random slope.

This is always true, in a way, with an interaction, but outside of LMMs we do not have random effects and so when we allow slopes to differ for different groups people, we do not know what an individual person’s slope was without the interaction so have no reference point. Put differently, outside of LMMs, in regular GLMs, we always assume that \(u_{2j} = 0\) so we wouldn’t really think about ‘predicting’ it when it is fixed at 0. In GLMs we only focus on how we allow the average slope to differ by level of the moderator. In LMMs we can interpret it the same way or we can interpret the moderator as ‘predicting’ the random slope.

3 Continuous Interactions in R

Aside from the notes about some minor interpretation differences, in general interactions in LMMs are analysed, graphed, and interpretted the same way as for GLMs.

First to avoid any issues around diagnostics etc. from haven labeled type data, we will convert the variable we are going to work with to numeric. Then we fit a LMM with an interaction between stress and neuroticism, mood as the outcome and a random intercept as the only random effect.

dm[, mood := as.numeric(mood)]
dm[, stress := as.numeric(stress)]
dm[, neuroticism := as.numeric(neuroticism)]

m <- lmer(mood ~ neuroticism * stress + (1 | ID), data = dm)

A quick check of the model diagnostics suggests that there may be an outlier on the residuals, but it is not too extreme and otherwise the data look fairly good. The residuals do not appear to follow a normal distribution that closely, partly due to the long left tail.

plot(modelDiagnostics(m), nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

Applying transformations to left skewed data is more difficult as generally transformations work on long right tails. A solution is to reverse the variable, apply the transformation and then again reverse it so that the direction is the same as it originally was. We could try a square root transformations which is milder than a log transformation. To reverse it, we subtract the variable from the sum of its minimum and maximum. Next we take its square root, then we reverse by again subtracting from the sum of its minimum and maximum, but square root transformed.

max(dm$mood) + min(dm$mood)
## [1] 8
## transform
dm[, moodtrans := sqrt(8) - sqrt(8 - mood)]

m <- lmer(moodtrans ~ neuroticism * stress + (1 | ID), data = dm)

md <- modelDiagnostics(m) 
plot(md, nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

The transformation appears to have modestly helped the distribution of residuals. Its not that clear whether it was bad enough to begin with and whether the transformation improved it enough that it is worth the difficulty in interpretation (mood is now square root transformed and that must be incorporated into its interpretation). For the lecture, we will proceed with the transformation, but in practice, consider whether this is worth it or only adds difficulty in understanding without improving / changing results much. We still see an extreme value present. We can remove that, as we have discussed in depth in previous lectures, update the model, and re-run diagnostics.

dmnoev <- dm[-md$extremeValues$Index]
mnoev <- update(m, data = dmnoev)
md <- modelDiagnostics(mnoev)
plot(md, nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

After removing one extreme value, another point is classified as an extreme value. We could remove that too. Note that its important that we remove it from the dmnoev dataset, not the original dm dataset, because the row numbers from these two datasets do not match and we are only removing the additional “new” extreme value found. It is relatively rare to actually iterate through this process of identifying and removing extreme values multiple times, but I do it to show how it woudl work if you wanted to.

After removing the additional extreme value, we would update the model again, and check diagnostics again.

dmnoev <- dmnoev[-md$extremeValues$Index]
mnoev <- update(mnoev, data = dmnoev)
md <- modelDiagnostics(mnoev)
plot(md, nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

Finally no extreme values are identified and the distributions are approximately normal. In practice it could take additional rounds of extreme value removal or you may decide to stop at one round.

At this point, we have a “clean” model and we can look at a summary of it. Although we have used summary() a lot in the past, we’ll introduce another function to help look at lmer() model results, modelTest(). In this lecture, we will only learn and interpret part of its output, with the rest of the output from modelTest() covered later. In addition to get nicely formatted results rather than a set of datasets containing the results, we use the APAStyler() function.

APAStyler(modelTest(mnoev))
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
##                           Term                    Est           Type
##  1:                (Intercept)  1.03*** [ 0.54, 1.53]  Fixed Effects
##  2:                neuroticism     0.06 [-0.01, 0.12]  Fixed Effects
##  3:         neuroticism:stress   -0.02* [-0.04, 0.00]  Fixed Effects
##  4:                     stress     0.08 [-0.05, 0.22]  Fixed Effects
##  5:          sd_(Intercept)|ID                   0.19 Random Effects
##  6:                      sigma                   0.26 Random Effects
##  7:                   Model DF                      6  Overall Model
##  8:                 N (Groups)                ID (50)  Overall Model
##  9:           N (Observations)                    181  Overall Model
## 10:                     logLik                 -37.53  Overall Model
## 11:                        AIC                  87.05  Overall Model
## 12:                        BIC                 106.24  Overall Model
## 13:                Marginal R2                   0.17  Overall Model
## 14:                Marginal F2                   0.21  Overall Model
## 15:             Conditional R2                   0.45  Overall Model
## 16:             Conditional F2                   0.83  Overall Model
## 17:        neuroticism (Fixed)   0.00/-0.01, p = .110   Effect Sizes
## 18:             stress (Fixed)   0.01/ 0.01, p = .231   Effect Sizes
## 19: neuroticism:stress (Fixed)   0.03/ 0.03, p = .015   Effect Sizes

The results show the regression coefficients, asterisks for p-values, and 95% confidence intervals in brackets for the fixed effects, the standard deviations of the random effects, the model degrees of freedom, which is how many parameters were estimated in the model total, and the number of people and observations. For now, we will ignore all the output under row 9, N (Observations). In the case of this model we can see the following.

A LMM was fit with 181 observations from 50 people. There was a significant neuroticism x stress interaction, b [95% CI] = -0.02 [-.04, .00], p < .05.

We can also pass multiple model results in a list together, which puts the results side by side. This is particularly helpful for comparing models with and without covariates or to evaluate whether removing extreme values changed the results substantially.

APAStyler(list(
  EV = modelTest(m),
  NOEV = modelTest(mnoev) ))
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
##                           Term                     EV                   NOEV
##  1:                (Intercept)  1.03*** [ 0.52, 1.54]  1.03*** [ 0.54, 1.53]
##  2:                neuroticism     0.06 [-0.02, 0.13]     0.06 [-0.01, 0.12]
##  3:         neuroticism:stress   -0.02* [-0.04, 0.00]   -0.02* [-0.04, 0.00]
##  4:                     stress     0.08 [-0.06, 0.22]     0.08 [-0.05, 0.22]
##  5:          sd_(Intercept)|ID                   0.19                   0.19
##  6:                      sigma                   0.28                   0.26
##  7:                   Model DF                      6                      6
##  8:                 N (Groups)                ID (50)                ID (50)
##  9:           N (Observations)                    183                    181
## 10:                     logLik                 -46.99                 -37.53
## 11:                        AIC                 105.97                  87.05
## 12:                        BIC                 125.23                 106.24
## 13:                Marginal R2                   0.15                   0.17
## 14:                Marginal F2                   0.18                   0.21
## 15:             Conditional R2                   0.41                   0.45
## 16:             Conditional F2                   0.69                   0.83
## 17:        neuroticism (Fixed)    0.00/0.00, p = .128   0.00/-0.01, p = .110
## 18:             stress (Fixed)    0.01/0.01, p = .266   0.01/ 0.01, p = .231
## 19: neuroticism:stress (Fixed)    0.03/0.03, p = .024   0.03/ 0.03, p = .015
##               Type
##  1:  Fixed Effects
##  2:  Fixed Effects
##  3:  Fixed Effects
##  4:  Fixed Effects
##  5: Random Effects
##  6: Random Effects
##  7:  Overall Model
##  8:  Overall Model
##  9:  Overall Model
## 10:  Overall Model
## 11:  Overall Model
## 12:  Overall Model
## 13:  Overall Model
## 14:  Overall Model
## 15:  Overall Model
## 16:  Overall Model
## 17:   Effect Sizes
## 18:   Effect Sizes
## 19:   Effect Sizes

These results show us that we have effectively the same results with or without those two extreme values we omitted. Indeed the only changes really observed are small changes in the 95% confidence intervals and a slight decrease in the residual standard deviation, sigma, with the extreme values removed (outside, of course of the models being based on 181 instead of 183 people). In a case like this, it will not make any real difference to the interpretation if the model with extreme values or without extreme values is used.

We can understand what the interaction means in the same way as for GLMs. The following 3D figure shows the slight warp in the regression plane between stress and neuroticism with (transformed) mood. The warping is because of the interaction.

## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

3.1 Plotting

Typically, 3D graphs are not used to plot interactions. Instead, a few exemplar lines are graphed showing the slope of one variable with the outcome at different values of the moderator. As with GLMs, we can use the visreg() function. Here, we’ll use neuroticism as the moderator. A common approach to picking level of the moderator is to use the Mean - 1 SD and Mean + 1 SD. To do that, we first need the mean and standard deviation of neuroticism, which we can get using egltable() after excluding duplicates by ID, since neuroticism only varies between units.

egltable("neuroticism", data = dmnoev[!duplicated(ID)])
##                     M (SD)
## 1: neuroticism 7.20 (2.08)
visreg(mnoev, xvar = "stress",
       by = "neuroticism", overlay=TRUE,
       breaks = c(7.2 - 2.08, 7.20 + 2.08),
       partial = FALSE, rug = FALSE)

The results show, in an easier to interpret way, what the negative interaction coefficient of \(b = -.02\) means, people with higher levels of neuroticism are more sensitive to the effects of stress. People lower in neuroticism are relatively less sensitive to the effects of stress, although in both cases, higher stress is associated with lower (transformed) mood.

Another common way of picking some exemplar values is to use the 25th and 75th percentiles. These work particularly well for very skewed distributions where the mean +/- SD could be outside the observed range of the data. Again we exclude duplicates by ID and then use the quantile() function to get the values, 6, and 9 for the 25th and 75th percentiles.

quantile(dmnoev[!duplicated(ID), neuroticism], na.rm = TRUE)
##   0%  25%  50%  75% 100% 
##    2    6    8    9   10
visreg(mnoev, xvar = "stress",
       by = "neuroticism", overlay=TRUE,
       breaks = c(6, 9),
       partial = FALSE, rug = FALSE)

3.2 Simple Effects

When working with models that have interactions, a common aid to interpretation is to test the simple effects / slopes from the model. For example, previously we graphed the association between stress and transformed mood at M - 1 SD and M + 1 SD on neuroticism. However, although visually both lines appeared to have a negative slope, we do not know from the graph alone whether there is a significant association betwee stress and transformed mood at both the low (M - 1 SD) and high (M + 1 SD) levels of neuroticism. To answer that, we need to test the simple slope of stress at specific values of neuroticism. Our default model does actually give us one simple slope: it is the simple slope of stress when \(neuroticism = 0\). However, as we can tell from the mean and standard deviation of neuroticism, 0 is very far outside the plausible range of values so that simple slope given to us by default from the model is not too useful. We could either center neuroticism and re-run the model, which would get us a different simple slope, or use post hoc functions to calculate simple slopes.

We will use the emtrends() function from the emmeans package to test the simple slopes. This function also works with GLMs, for your reference.

The emtrends() function take a model as its first argument, then the variable that you want to calculate a simple slope for, here stress, the argument at requires a list of specific values of the moderator, and then we tell it how we want degrees of freedom calculated (note this only applies to lmer models). We store the results in an R object, mem and then call summary() to get a summary table. The infer = TRUE argument is needed to summary() if you want p-values.

mem <- emtrends(mnoev, var = "stress",
                at = list(neuroticism = c(7.2 - 2.08, 7.2 + 2.08)),
                lmer.df = "satterthwaite")

summary(mem, infer=TRUE)
##  neuroticism stress stress.trend     SE  df lower.CL upper.CL t.ratio p.value
##         5.12   3.82      -0.0311 0.0273 171   -0.085   0.0229 -1.137  0.2572 
##         9.28   3.82      -0.1237 0.0252 177   -0.173  -0.0739 -4.903  <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

The relevant parts of the output, for us, are the columsn for stress.trend which are the simple slopes, the values of neuroticism which tell us at what values of neuroticism we have calculated simple slopes, the confidence intervals, lower.CL and upper.CL, 95% by default, and the p-value. From these results, we can see that when \(neuroticism = 5.12\) there is no significant association between stress and transformed mood, but there is when \(neuroticism = 9.28\).

3.3 Sample Write Up

With all of this information, we can plan out some final steps for a polished write up of the results. First, let’s get exact p-values for all our results. We can do this through options to pcontrol in APAStyler(). We also re-print the simple slopes here.

APAStyler(list(
  EV = modelTest(m),
  NOEV = modelTest(mnoev)),
  pcontrol = list(digits = 3, stars = FALSE, includeP = TRUE,
                  includeSign = TRUE, dropLeadingZero = TRUE))
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
##                           Term                          EV
##  1:                (Intercept)  1.03p < .001 [ 0.52, 1.54]
##  2:                neuroticism  0.06p = .133 [-0.02, 0.13]
##  3:         neuroticism:stress -0.02p = .026 [-0.04, 0.00]
##  4:                     stress  0.08p = .273 [-0.06, 0.22]
##  5:          sd_(Intercept)|ID                        0.19
##  6:                      sigma                        0.28
##  7:                   Model DF                           6
##  8:                 N (Groups)                     ID (50)
##  9:           N (Observations)                         183
## 10:                     logLik                      -46.99
## 11:                        AIC                      105.97
## 12:                        BIC                      125.23
## 13:                Marginal R2                        0.15
## 14:                Marginal F2                        0.18
## 15:             Conditional R2                        0.41
## 16:             Conditional F2                        0.69
## 17:        neuroticism (Fixed)         0.00/0.00, p = .128
## 18:             stress (Fixed)         0.01/0.01, p = .266
## 19: neuroticism:stress (Fixed)         0.03/0.03, p = .024
##                            NOEV           Type
##  1:  1.03p < .001 [ 0.54, 1.53]  Fixed Effects
##  2:  0.06p = .115 [-0.01, 0.12]  Fixed Effects
##  3: -0.02p = .017 [-0.04, 0.00]  Fixed Effects
##  4:  0.08p = .238 [-0.05, 0.22]  Fixed Effects
##  5:                        0.19 Random Effects
##  6:                        0.26 Random Effects
##  7:                           6  Overall Model
##  8:                     ID (50)  Overall Model
##  9:                         181  Overall Model
## 10:                      -37.53  Overall Model
## 11:                       87.05  Overall Model
## 12:                      106.24  Overall Model
## 13:                        0.17  Overall Model
## 14:                        0.21  Overall Model
## 15:                        0.45  Overall Model
## 16:                        0.83  Overall Model
## 17:        0.00/-0.01, p = .110   Effect Sizes
## 18:        0.01/ 0.01, p = .231   Effect Sizes
## 19:        0.03/ 0.03, p = .015   Effect Sizes
summary(mem, infer=TRUE)
##  neuroticism stress stress.trend     SE  df lower.CL upper.CL t.ratio p.value
##         5.12   3.82      -0.0311 0.0273 171   -0.085   0.0229 -1.137  0.2572 
##         9.28   3.82      -0.1237 0.0252 177   -0.173  -0.0739 -4.903  <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

Now we will make a polished, finalized figure. I have customized the colours, and turned off the legends. In place of legends, I have manually added text annotations including the simple slopes and confidence intervals and p-values for the simple slopes1 For your reference, it took about 6 trial and errors of different x and y values and angles to get the text to line up about right. I did not magically get the values to use to get a graph that I thought looked nice. That is why I think sometimes it is easier to add this sort of text after the fact in your slides of papers rather than building it into the code..

visreg(m, xvar = "stress",
       by = "neuroticism", overlay=TRUE,
       breaks = c(7.2 - 2.08, 7.20 + 2.08),
       partial = FALSE, rug = FALSE, gg=TRUE,
       xlab = "Daily Stress",
       ylab = "Predicted Daily Mood") +
  scale_color_manual(values = c("5.12" = "black", "9.28" = "grey70")) +
  theme_pubr() +
  guides(colour = FALSE, fill = FALSE) +
  annotate(geom = "text", x = 4.2, y = 1.21, label = "High Neuroticism: b = -.03 [-.09, .02], p = .26",
           angle = -10) + 
  annotate(geom = "text", x = 4.4, y = 1.045, label = "Low Neuroticism: b = -.12 [-.17, -.07], p < .001",
           angle = -33.8)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

A linear mixed model using restricted maximum likelihood was used to test whether the association of daily stress on daily mood is moderated by baseline neuroticism scores. All predictors were included as fixed effects and a random intercept by participant was included. Visual diagnostics showed that mood was not normally distributed, so it was square root transformed, referred to as transformed mood hereafter. Additionally, two outliers were identified on the residuals and these were excluded leaving a total of 181 observations from 50 people.

There was a significant daily stress x neuroticism interaction, such that higher neuroticism scores were associated with a more negative association between daily stress and same day transformed mood, b [95% CI] = -0.02 [-0.04, 0.00], p = .017. To help interpret the interaction, the simple slopes of daily stress on transformed mood were tested and graphed at low (M - 1 SD) and high (M + 1 SD) values of neuroticism (see Figure XX). The results revealed that there was a significant, negative association between daily stress and same day transformed mood at high, but not low levels of neuroticism.

Sensitivity analyses not excluding outliers revealed that results remained unchanged.

4 Continuous x Categorical Interactions in R

Continuous x Categorical interactions are conducted much as continuous x continuous interactions. Typically with continuous x categorical interactions, simple slopes for the continuous variable are calculated at all levels of the categorical variable.

Here we will work with a three level, dummy coded variable, stress with low being the reference and two dummy codes, one for mid and one for high stress. The outcome is transformed mood we used earlier.

## create a "categorical" stress variable
dm[, stress3 := cut(stress, breaks = quantile(stress, probs = c(0, 1/3, 2/3, 1), na.rm=TRUE),
                    labels = c("Low", "Mid", "High"),
                    include.lowest = TRUE)]

m <- lmer(moodtrans ~ neuroticism * stress3 + (1 | ID), data = dm)

The model diagnostics look relatively good, albeit not perfect.

plot(modelDiagnostics(m), nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

With reasonable diagnostics, we can look at a summary.

We have a “clean” model and we can look at a summary of it.

summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: moodtrans ~ neuroticism * stress3 + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 129.7
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.781 -0.617  0.139  0.530  2.767 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.0360   0.190   
##  Residual             0.0824   0.287   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)               1.2710     0.1619  85.4891    7.85  1.1e-11 ***
## neuroticism              -0.0103     0.0227  94.7196   -0.45     0.65    
## stress3Mid                0.1191     0.2141 175.4900    0.56     0.58    
## stress3High               0.3798     0.4819 165.3312    0.79     0.43    
## neuroticism:stress3Mid   -0.0289     0.0285 176.2966   -1.02     0.31    
## neuroticism:stress3High  -0.0830     0.0578 168.1476   -1.44     0.15    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) nrtcsm strs3M strs3H nrt:3M
## neuroticism -0.956                            
## stress3Mid  -0.541  0.525                     
## stress3High -0.230  0.234  0.243              
## nrtcsm:st3M  0.550 -0.585 -0.966 -0.255       
## nrtcsm:st3H  0.263 -0.291 -0.274 -0.983  0.308

Factor variables in interactions do not work currently with modelTest(), so if we wanted to use it, we’d need to manually dummy code the categorical variable. The results are identical.

dm[, stressHigh := as.integer(stress3 == "High")]
dm[, stressMid := as.integer(stress3 == "Mid")]

malt <- lmer(moodtrans ~ neuroticism * (stressMid + stressHigh) + (1 | ID), data = dm)

APAStyler(modelTest(malt))
## Parameters and CIs are based on REML, 
## but modelTests requires ML not REML fit for comparisons, 
## and these are used in effect sizes. Refitting.
##                               Term                    Est           Type
##  1:                    (Intercept)  1.27*** [ 0.95, 1.59]  Fixed Effects
##  2:                    neuroticism    -0.01 [-0.05, 0.03]  Fixed Effects
##  3:         neuroticism:stressHigh    -0.08 [-0.20, 0.03]  Fixed Effects
##  4:          neuroticism:stressMid    -0.03 [-0.08, 0.03]  Fixed Effects
##  5:                     stressHigh     0.38 [-0.56, 1.32]  Fixed Effects
##  6:                      stressMid     0.12 [-0.30, 0.54]  Fixed Effects
##  7:              sd_(Intercept)|ID                   0.19 Random Effects
##  8:                          sigma                   0.29 Random Effects
##  9:                       Model DF                      8  Overall Model
## 10:                     N (Groups)                ID (50)  Overall Model
## 11:               N (Observations)                    183  Overall Model
## 12:                         logLik                 -51.01  Overall Model
## 13:                            AIC                 118.02  Overall Model
## 14:                            BIC                 143.70  Overall Model
## 15:                    Marginal R2                   0.11  Overall Model
## 16:                    Marginal F2                   0.12  Overall Model
## 17:                 Conditional R2                   0.37  Overall Model
## 18:                 Conditional F2                   0.59  Overall Model
## 19:            neuroticism (Fixed)    0.01/0.00, p = .638   Effect Sizes
## 20:              stressMid (Fixed)    0.00/0.01, p = .578   Effect Sizes
## 21:             stressHigh (Fixed)    0.01/0.00, p = .419   Effect Sizes
## 22:  neuroticism:stressMid (Fixed)    0.01/0.01, p = .306   Effect Sizes
## 23: neuroticism:stressHigh (Fixed)    0.02/0.01, p = .144   Effect Sizes
##                               Term                    Est           Type

4.1 Plotting

These interactions are not significant, so often we would simply drop the interaction term and report main effects. For the sake of example only, we will plot anyways.

With continuous x categorical interactions, the easiest approach is to plot the simple slope of the continuous variable by the categorical one as shown in the following.

visreg(m, xvar = "neuroticism",
       by = "stress3", overlay=TRUE,
       partial = FALSE, rug = FALSE)

4.2 Simple Effects

These interactions are not significant, so often we would simply drop the interaction term and report main effects. For the sake of example only, we will continue to look at simple effects anyways.

When working with models that have interactions, a common aid to interpretation is to test the simple effects / slopes from the model. For example, previously we graphed the association between neuroticism and transformed mood at each level of the categorical stress3 variable. However, we cannot tell from the graph whether neuroticism is significantly associated with transformed mood for any specific level of stress. To answer that, we need to test the simple slope of neuroticism at specific values of neuroticism. Our default model does actually give us one simple slope: it is the simple slope of neuroticism when stress3 = low$, but we might want more.

We will use the emtrends() function from the emmeans package to test the simple slopes.

The emtrends() function take a model as its first argument, then the variable that you want to calculate a simple slope for, here stress, the argument at requires a list of specific values of the moderator, and then we tell it how we want degrees of freedom calculated (note this only applies to lmer models). We store the results in an R object, mem and then call summary() to get a summary table. The infer = TRUE argument is needed to summary() if you want p-values.

mem <- emtrends(m, var = "neuroticism",
                at = list(stress3 = c("Low", "Mid", "High")),
                lmer.df = "satterthwaite")

summary(mem, infer=TRUE)
##  neuroticism stress3 neuroticism.trend     SE    df lower.CL upper.CL t.ratio
##         7.34 Low               -0.0103 0.0227  94.7  -0.0552  0.03471 -0.453 
##         7.34 Mid               -0.0392 0.0239  95.7  -0.0866  0.00815 -1.643 
##         7.34 High              -0.0933 0.0556 175.5  -0.2031  0.01653 -1.677 
##  p.value
##  0.6514 
##  0.1036 
##  0.0954 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

The relevant parts of the output, for us, are the columsn for neuroticism.trend which are the simple slopes, the values of stress3 which tell us at what values of stress3 we have calculated simple slopes, the confidence intervals, lower.CL and upper.CL, 95% by default, and the p-value. From these results, we can see that neuroticism is not significantly associated with transformed mood for any stress category, although its closest for high stress.

5 Categorical x Categorical Interactions in R

Categorical x Categorical interactions are conducted comparably, although more contrasts / simple effect follow-ups are possible.

Here we will work with a three level, dummy coded variable, stress with low being the reference and two dummy codes, one for mid and one for high stress. The outcome is transformed mood we used earlier. We also work with a three level neuroticism variable.

## create a categorical neuroticism variable
dm[, neuro3 := cut(neuroticism, breaks = quantile(neuroticism, probs = c(0, 1/3, 2/3, 1), na.rm=TRUE),
                    labels = c("Low", "Mid", "High"),
                    include.lowest = TRUE)]

m <- lmer(moodtrans ~ neuro3 * stress3 + (1 | ID), data = dm)

The model diagnostics look relatively good, albeit not perfect, with only one modest extreme value.

plot(modelDiagnostics(m), nrow = 2, ncol = 2, ask = FALSE)
## `geom_smooth()` using formula 'y ~ x'

With reasonable diagnostics, we can look at a summary.

summary(m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: moodtrans ~ neuro3 * stress3 + (1 | ID)
##    Data: dm
## 
## REML criterion at convergence: 124.2
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -2.674 -0.657  0.147  0.526  2.807 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.0380   0.195   
##  Residual             0.0811   0.285   
## Number of obs: 183, groups:  ID, 50
## 
## Fixed effects:
##                        Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)              1.2033     0.0650  81.6731   18.51   <2e-16 ***
## neuro3Mid                0.0310     0.1175  89.8943    0.26    0.793    
## neuro3High              -0.0566     0.1202 106.2919   -0.47    0.639    
## stress3Mid              -0.0746     0.0783 172.9876   -0.95    0.342    
## stress3High             -0.3247     0.1493 160.4072   -2.18    0.031 *  
## neuro3Mid:stress3Mid     0.0192     0.1325 173.9233    0.14    0.885    
## neuro3High:stress3Mid   -0.0936     0.1314 173.1645   -0.71    0.477    
## neuro3Mid:stress3High    0.3122     0.2269 158.1092    1.38    0.171    
## neuro3High:stress3High  -0.1278     0.2039 164.9928   -0.63    0.532    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ner3Md nr3Hgh strs3M strs3H n3M:3M n3H:3M n3M:3H
## neuro3Mid   -0.554                                                 
## neuro3High  -0.541  0.300                                          
## stress3Mid  -0.590  0.327  0.319                                   
## stress3High -0.270  0.149  0.146  0.279                            
## nr3Md:str3M  0.348 -0.640 -0.189 -0.591 -0.164                     
## nr3Hgh:st3M  0.351 -0.195 -0.676 -0.596 -0.166  0.352              
## nr3Md:str3H  0.178 -0.340 -0.096 -0.183 -0.658  0.359  0.109       
## nr3Hgh:st3H  0.198 -0.109 -0.433 -0.204 -0.732  0.120  0.444  0.482

5.1 Plotting

With categorical x categorical interactions, visreg() produces OK but not great figures as shown in the following. We can see the means of transformed mood for all 9 cells (the cross of 3 level of neuroticism x 3 levels of stress).

visreg(m, xvar = "neuro3",
       by = "stress3", overlay=TRUE,
       partial = FALSE, rug = FALSE)

5.2 Simple Effects

When working with two categorical interactions (or, as an aside, with a categorical predictor with >2 levels where you want to test various group differences), the emmeans() function from the emmeans package is helpful. We can get the means of each neuroticism group by stress group and get confidence intervals and p-values. These p-values test whether each mean is different from zero, by default.

em <- emmeans(m, "neuro3", by = "stress3",
              lmer.df = "satterthwaite")
summary(em, infer = TRUE)
## stress3 = Low:
##  neuro3 emmean     SE    df lower.CL upper.CL t.ratio p.value
##  Low     1.203 0.0650  81.7    1.074    1.333 18.505  <.0001 
##  Mid     1.234 0.0978  93.8    1.040    1.429 12.620  <.0001 
##  High    1.147 0.1011 118.0    0.947    1.347 11.345  <.0001 
## 
## stress3 = Mid:
##  neuro3 emmean     SE    df lower.CL upper.CL t.ratio p.value
##  Low     1.129 0.0659  84.4    0.998    1.260 17.119  <.0001 
##  Mid     1.179 0.0842  70.3    1.011    1.347 13.999  <.0001 
##  High    0.979 0.0775  71.7    0.824    1.133 12.626  <.0001 
## 
## stress3 = High:
##  neuro3 emmean     SE    df lower.CL upper.CL t.ratio p.value
##  Low     0.879 0.1458 173.4    0.591    1.167  6.024  <.0001 
##  Mid     1.222 0.1608 173.9    0.904    1.539  7.596  <.0001 
##  High    0.694 0.1164 152.3    0.464    0.924  5.962  <.0001 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

A nice plot, with confidence intervals for the fixed effects, can be obtained by using the emmip() function from the emmeans package. It takes as input the results from emmeans(), not the lmer() model results directly. Here is a simple plot showing the categorical interactions. Note that with this approach, you could basically fit the same model(s) that you would with a repeated measures or mixed effects ANOVA model, with the advantage that LMMs do not require balanced designs and allow both categorical and continuous predictors (e.g., you could include continuous covariates easily). GLMs and (G)LMMs can do everything that t-tests and various ANOVAs can, but with greater flexibility.

emmip(em, stress3 ~ neuro3, CIs = TRUE) +
  theme_pubr() +
  ylab("Predicted Transformed Mood")

If you want pairwise comparisons, you can get all possible pairwise comparisons between neuroticism levels by stress or stress levels by neuroticism using the pairs() function.

## pairwise comparisons of neuroticism by stress
pairs(em, by = "stress3")
## stress3 = Low:
##  contrast   estimate    SE    df t.ratio p.value
##  Low - Mid   -0.0310 0.117  89.9 -0.264  0.9625 
##  Low - High   0.0566 0.120 106.3  0.471  0.8851 
##  Mid - High   0.0875 0.141 105.8  0.622  0.8082 
## 
## stress3 = Mid:
##  contrast   estimate    SE    df t.ratio p.value
##  Low - Mid   -0.0501 0.107  75.3 -0.469  0.8861 
##  Low - High   0.1502 0.102  76.8  1.476  0.3079 
##  Mid - High   0.2004 0.114  70.9  1.751  0.1938 
## 
## stress3 = High:
##  contrast   estimate    SE    df t.ratio p.value
##  Low - Mid   -0.3432 0.217 173.9 -1.581  0.2568 
##  Low - High   0.1844 0.187 171.9  0.988  0.5853 
##  Mid - High   0.5276 0.199 170.6  2.657  0.0234 
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 3 estimates
## pairwise comparisons of stress by neuroticism
pairs(em, by = "neuro3")
## neuro3 = Low:
##  contrast   estimate     SE  df t.ratio p.value
##  Low - Mid    0.0746 0.0783 173  0.953  0.6073 
##  Low - High   0.3247 0.1493 160  2.175  0.0786 
##  Mid - High   0.2501 0.1480 157  1.690  0.2122 
## 
## neuro3 = Mid:
##  contrast   estimate     SE  df t.ratio p.value
##  Low - Mid    0.0554 0.1070 174  0.518  0.8626 
##  Low - High   0.0125 0.1709 156  0.073  0.9971 
##  Mid - High  -0.0429 0.1598 145 -0.269  0.9610 
## 
## neuro3 = High:
##  contrast   estimate     SE  df t.ratio p.value
##  Low - Mid    0.1682 0.1056 170  1.594  0.2510 
##  Low - High   0.4526 0.1389 169  3.257  0.0039 
##  Mid - High   0.2843 0.1146 150  2.480  0.0377 
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 3 estimates

You can also get custom contrasts. For example, if we wanted to compare high neuroticism to the average of low and mid neuroticism at each level of stress (H1) and secondly see if low and mid neuroticism differ (H2). The list gives the specific contrast weights, which are directly applied to the means we saw earlier from emmeans().

contrast(em,
         list(
           H1 = c(.5, .5, -1),
           H2 = c(1, -1, 0)),
         by = "stress3")
## stress3 = Low:
##  contrast estimate     SE    df t.ratio p.value
##  H1         0.0721 0.1169 110.4  0.617  0.5388 
##  H2        -0.0310 0.1174  89.9 -0.264  0.7927 
## 
## stress3 = Mid:
##  contrast estimate     SE    df t.ratio p.value
##  H1         0.1753 0.0942  72.8  1.862  0.0667 
##  H2        -0.0501 0.1070  75.3 -0.469  0.6405 
## 
## stress3 = High:
##  contrast estimate     SE    df t.ratio p.value
##  H1         0.3560 0.1592 167.7  2.236  0.0266 
##  H2        -0.3432 0.2171 173.9 -1.581  0.1158 
## 
## Degrees-of-freedom method: satterthwaite

To help you see directly how the contrast weights are applied, we can apply them directly to the estimated means.

## all the means
as.data.frame(em)
##   neuro3 stress3 emmean      SE     df lower.CL upper.CL
## 1    Low     Low 1.2033 0.06503  81.67   1.0740   1.3327
## 2    Mid     Low 1.2343 0.09781  93.78   1.0401   1.4285
## 3   High     Low 1.1467 0.10108 118.05   0.9466   1.3469
## 4    Low     Mid 1.1287 0.06593  84.39   0.9976   1.2598
## 5    Mid     Mid 1.1789 0.08421  70.26   1.0109   1.3468
## 6   High     Mid 0.9785 0.07750  71.72   0.8240   1.1330
## 7    Low    High 0.8786 0.14584 173.36   0.5908   1.1665
## 8    Mid    High 1.2218 0.16085 173.95   0.9043   1.5393
## 9   High    High 0.6942 0.11644 152.34   0.4641   0.9242
## just the first three, all for stress3 = low
as.data.frame(em)$emmean[1:3]
## [1] 1.203 1.234 1.147
## apply H1 weights
sum(as.data.frame(em)$emmean[1:3] * c(.5, .5, -1))
## [1] 0.07207
## apply H2 weight
sum(as.data.frame(em)$emmean[1:3] * c(1, -1, 0))
## [1] -0.03095

You could get even more specific hypotheses about group differences, if you wanted using by = NULL all the means are in a row of 9 and we can apply weights to each. For example, here we contrast the average of low neuroticism at low stress and low neuroticism at mid stress to high neuroticism and high stress (H1). For more examples of weighting schemes for many different possible specific contrasts, wee: https://stats.idre.ucla.edu/spss/faq/how-can-i-test-contrasts-and-interaction-contrasts-in-a-mixed-model/ .

## all the means
as.data.frame(em)
##   neuro3 stress3 emmean      SE     df lower.CL upper.CL
## 1    Low     Low 1.2033 0.06503  81.67   1.0740   1.3327
## 2    Mid     Low 1.2343 0.09781  93.78   1.0401   1.4285
## 3   High     Low 1.1467 0.10108 118.05   0.9466   1.3469
## 4    Low     Mid 1.1287 0.06593  84.39   0.9976   1.2598
## 5    Mid     Mid 1.1789 0.08421  70.26   1.0109   1.3468
## 6   High     Mid 0.9785 0.07750  71.72   0.8240   1.1330
## 7    Low    High 0.8786 0.14584 173.36   0.5908   1.1665
## 8    Mid    High 1.2218 0.16085 173.95   0.9043   1.5393
## 9   High    High 0.6942 0.11644 152.34   0.4641   0.9242
contrast(em,
         list(
           H1 = c(.5, 0, 0, .5, 0, 0, 0, 0, -1)),
         by = NULL)
##  contrast estimate    SE  df t.ratio p.value
##  H1          0.472 0.128 130 3.694   0.0003 
## 
## Degrees-of-freedom method: satterthwaite

6 Summary

6.1 Conceptual

Key points to take away conceptually are:

6.2 Code

Function What it does
lmer() estimate a LMM
confint() calculate confidence intervals for a LMM
visreg() create marginal or conditional graphs from a LMM
modelDiagnostics() evaluate model diagnostics for LMMs including of multivariate normality
summary() get a summary of model results
modelTest() along with APAStyler() get a nicely formatted summary of a model results.
emmeans() test specific means from a model.
emtrends() test simple slopes from a model.
contrast() test custom contrasts on a set of means from emmeans().